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Abstract. We describe a numerical method to simulate an elastic shell im- 
mersed in a viscous incompressible fluid. The method is developed as an 
extension of the immersed boundary method using shell equations based on 
the Kirchhoff-Love and the planar stress hypotheses. A detailed derivation of 
the shell equations used in the numerical method is presented. This derivation 
as well as the numerical method, use techniques of differential geometry in 
an essential way. Our main motivation for the development of this method 
is its use in the construction of a comprehensive three-dimensional computa- 
tional model of the cochlea (the inner ear). The central object of study within 
the cochlea is the "basilar membrane", which is immersed in fluid and whose 
elastic properties rather resemble those of a shell. We apply the method to a 
specific example, which is a prototype of a piece of the basilar membrane and 
study the convergence of the method in this case. Some typical features of 
cochlear mechanics are already captured in this simple model. In particular, 
numerical experiments have shown a traveling wave propagating from the base 
to the apex of the model shell in response to external excitation in the fluid. 



1. Introduction 

This paper describes a general method of simulation of an elastic shell immersed 
in a viscous incompressible fluid. The method is developed as an extension of 
the immersed boundary method originally introduced by Peskin and McQueen to 
study the blood flow in the heart. The immersed boundary method has proved to 
be particularly useful for computer simulation of various biofluid dynamic systems. 
In this framework the elastic (and possibly active) biological tissue is treated as a 
collection of elastic fibers immersed in a viscous incompressible fluid. This formu- 
lation of the method together with references to many applications can be found in 
. A partial list of the applications of the immersed boundary method includes 
in addition to the blood flow in the heart (see the extensive work of Peskin and 
McQueen, e. g. ^JEU) also platelet aggregation during blood clotting 0|, flow 
of suspensions |SJ I15j , aquatic animal locomotion 0] , a two-dimensional model of 
cochlear fluid mechanics Q and flow in collapsible tubes [HJ. For a recent review of 
the immersed boundary method, see |12j . 

Many man-made materials can be modeled as elastic shells and elastic shells 
are also ubiquitous in nature, however our motivation for the development of the 
method comes from the study of the cochlea. The auditory signal processing in 
the cochlea depends crucially on the dynamics of the basilar membrane, which is 
immersed in a viscous incompressible fluid of the cochlea, and despite its name the 
basilar membrane is actually an elastic shell. The numerical method for elastic 
shell-fluid interaction presented here was subsequently used in the construction of 

Computation was performed at the Pittsburgh Supercomputing Center under allocation 
MCA93S004 from the MetaCenter Allocations Committee. 
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a complete three-dimensional computational model of the macro-mechanics of the 
cochlea which incorporates the intricate curved cochlear anatomy. The results of 
this work will be reported in future publications. 

The cochlea is the part of the inner ear where sound waves are transformed into 
electrical pulses which are carried by neurons to the brain. It is a small snail- 
shell-like cavity in the temporal bone, which has two openings, the oval window 
and the round window. The cavity is filled with fluid, which is sealed in by two 
elastic membranes covering these windows. The spiral canal of the cochlea is divided 
lengthwise by the long and narrow basilar membrane into two passages that connect 
with each other at the apex. External sounds set the ear drum in motion, which 
is conveyed to the inner ear by three small bones of the middle ear. These bones 
function as an impedance matching device, focusing the energy of the ear drum on 
the oval window of the cochlea. This piston-like motion against the oval window 
displaces the fluid of the cochlea generating traveling waves that propagate along 
the basilar membrane. The vibrations of the basilar membrane are detected by 
thousands of microscopic sensory receptors, called hair cells, located on the surface 
of the basilar membrane. The auditory signal processing in the cochlea is completed 
by the hair cells converting these mechanical stimuli into action potentials in the 
neurons attached to them, relaying this information to the brain. 

Practically everything we know about the passive wave propagation in the cochlea 
was discovered in the 1940s by Georg von Bekesy who carried out experiments in 
cochleae extracted from human cadavers. Von Bekesy observed that a pure tone 
input sound generates a traveling wave which reaches its peak at a speciffic loca- 
tion along the basilar membrane exciting only a narrow band of hair cells. This 
characteristic location depends on the tone's frequency. By pinching the basilar 
membrane with a tiny probe and observing the resulting displacement von Bekesy 
discovered that the basilar membrane is, in fact, not a membrane, i. e. it is not 
under inner tension, but an elastic shell, whose compliance varies exponentially 
along the membrane. Von Bekesys extensive experimental work is summarized in 
his book "Experiments in Hearing" For an excellent summary of more recent 
work on the cochlea, see UJ. 

Cochlear mechanics has been an active area of research ever since von Bekesy's 
fundamental contributions, yet many important questions are still open. Presently 
there is no complete understanding of the mechanisms responsible for the extreme 
sensitivity, sharp frequency selectivity and broad dynamic range of hearing. The 
most rigorous mathematical analysis of the cochlea was carried out by Leveque, Pe- 
skin and Lax in 9 . In their model the cochlea is represented by a two-dimensional 
plane (i.e. a strip of infinite length and infinite depth) and the basilar membrane, by 
a straight line of harmonic oscillators dividing the fluid plane into two halves. The 
linearized equations are reduced to a functional equation by applying the Fourier 
transform in the direction parallel to the basilar membrane and then solving the 
resulting ordinary differential equations in the normal direction. The functional 
equation derived in this way is solved analytically, and the solution is evaluated 
both numerically and also asymptotically (by the method of stationary phase). 
This analysis reveals that the waves in the cochlea resemble shallow water waves, 
i.e. ripples on the surface of a pond. A distinctive feature of this paper is the 
(then speculative) consideration of negative basilar membrane friction, i.e., of an 
amplification mechanism operating within the cochlea. 
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Like the cochlea, even the simplest three-dimensional fluid-shell systems appear 
to be too difficult to analyze, but using the methods presented here it is possible to 
construct computational models for such systems. The construction of computa- 
tional models for fluid-shell configurations is however not easy, and the development 
of these methods required supercomputing resources. Nevertheless, with the rapid 
advances in computer technology, such computations will soon be feasible on a 
workstation. The cochlea example is important also because the fluid-shell system 
simulation in this case can be compared with the extensive body of theoretical and 
experimental research. 

The theory of elastic plates and shells is a classical mathematical subject, which 
is also an active area of contemporary research (see for example |21 1111 ITTj ). For 
completeness, we derive in section the elastic shell equations upon which the 
numerical method is based. Shell theory is very naturally described in the language 
of differential geometry and it is perhaps in this respect that the presentation in 
section |21 somewhat differs from other accounts of the subject. Imagine a material 
composed of very rigid straight line segments (fibers) coupled together and assume 
that these fibers are perpendicular to some imaginary surface S in the middle of 
the material. In other words, the material has a structure of a normal vector 
bundle over the surface S. We shall assume that the base surface S is free to 
undergo arbitrary small elastic deformations. The Kirchhoff-Love hypothesis, is 
the assumption that the deformation of the bundle is such that the fibers of the 
deformed bundle are perpendicular to the deformed middle surface, and that these 
fibers are not stretched during the deformation. In the language of differential 
geometry this means that the shell has the structure of a normal vector bundle 
which is preserved under elastic deformations. We shall call such a material an 
elastic bundle. The Kirchhoff-Love hypothesis implies that a deformation of the 
base surface completely determines the deformation of the whole bundle. To ensure 
that the fibers do not intersect each other we must assume that their length, i.e., 
the thickness of the bundle, is smaller than the radius of maximal curvature of the 
base surface. 

Taking the three-dimensional linear theory of elasticity as our starting point, and 
using the Kirchhoff-Love hypothesis, it is possible to derive the equations which 
completely determine small deformations of elastic bundles. It turns out how- 
ever that a realistic model has to satisfy the hypothesis of planar stress. Strictly 
speaking, in linear elasticity the planar stress hypothesis is not consistent with the 
Kirchhoff-Love assumption. We utilize a common approach to "reconciling" the 
two assumptions (e.g see 7 ). The resulting equations constitute a system of three 
partial differential equations in three unknown functions. They express the elastic 
force exerted by the shell as a linear fourth order differential operator applied to the 
displacement vector field. This differential operator is intrinsic to the base surface 
of the bundle and its coefficients are tensorial quantities determined by the geom- 
etry and the elastic properties of the bundle. This intrinsic geometric formulation 
is essential in the formulation of the numerical method described in the following 
section. 

Section |21 outlines the immersed boundary method for shells. This is a modifi- 
cation of the immersed boundary method as described in |13j . The main difference 
in the algorithm is in the computation of the force that the material applies to the 
fluid. Here it is computed by discretizing the shell equations described in section^ 
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We describe in detail the discrete differential operators defined on the shell surface 
which are used in the force computation. 

In section 0] we describe a test model in which the shell resembles a piece of a 
basilar membrane immersed in fluid. We study the convergence of the algorithm in 
this case and describe the numerical experiments carried out with this model. These 
experiments have reproduced some of the typical features of cochlear mechanics, 
such as the traveling wave propagating along the basilar membrane in response to 
external excitation of the fluid. 

This work is based on the author's Ph.D. thesis completed at the Courant In- 
stitute (NYU) under the supervision of Professor Charles S. Peskin. I would like 
to thank Professor Peskin for his support, encouragement and patient guidance. I 
would also like to thank Professor David McQueen for many conversations in which 
I learnt a lot about scientific computing. I'd like to thank Professor Karl Grosh for 
helping me understand shell theory. ' 

2. Linear Elastic Deformations of Vector Bundles (Shell Theory) 

Let M C R 3 be a normal bundle over a surface S. We describe S by a coordinate 
chart <f> : f2 — > R 3 , where il is a domain in R 2 . Let n : Cl — » M be the unit normal 
vector field on S. The natural chart for M is <£> : ft x (—ho,ho)^M given by 

®(qi,q2,t) = <t>{qi,q2) +tn(q 1 ,q 2 ), 

and we assume, for simplicity, that the thickness of the bundle (= 2ho) is constant. 
Let $ e : x (—ho, ho) — > R 3 describe a 1-parameter family of deformations of 
M, such that <E>o = Our basic assumption is that each <E> e preserves the bundle 
structure of M . Clearly, such a deformation is completely determined by the de- 
formation of the base space S. In the framework of linear elasticity we work with 
infinitesimal deformations, i.e., vector fields. We will assume that S is free to un- 
dergo an arbitrary infinitesimal deformation ip — ^ | e= o0e and we shall determine 
the corresponding infinitesimal deformation V — ^| e= o$ e of M. 

We begin with some preliminary remarks about the geometry of M . Throughout 
this chapter the indices are raised and lowered with respect to h, the metric of S. 
We adopt the convention that greek indices take the values 1,2, while roman indices 
run through 1, 2, 3. The standard metric on R 3 will be denoted by 6: 

S(X, Y) =< X, Y >= X ■ Y, 

d a = -^j- is a partial derivative and V is the standard flat connection of R 3 . 

2.1. The Geometry of the Vector Bundle. Let h = <\>* (&\s) denote the metric 
of S in the chart <fi. The components of h are 

h a/ 3 = d a (f> ■ dpcf). 

Throughout this paper greek tensor indices are raised and lowered with respect to 
the metric h. Let N p denote the unit normal vector to S at the point p G M. 
Thus, n(<7i,<72) = ^<f>(qi.q 2 )- We shall often abuse notation and identify n with 
N. Similarly for other vector fields. The second fundamental form of S is the 
symmetric bilinear form acting on vectors tangent to the surface S defined by 

b(X, Y) =< V X N, Y > . 
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Its components in the chart (j> are 

b a fi = d a n ■ dfi4>. 

The central role in the following analysis is played by the map 9 = 9 t , which we 
define to be the flow of N: 

9t(p)=P + tN p , pes. 

The differential of 9 is the map 

: T p S — > Tg(p)St, 

where St — 9 t (S). Let X G T p S and let a be a curve on S such that <r(0) = p, 
cr'(O) = X. Then 

9*X = 9 ^ a 1 

= (9oa)' 

= (a + tNoa)' 

= X + tVxN 

We shall regard 9* as a symmetric tensor field on S. The components 9 a p of 9* can 
be expressed in the chart <f> as follows: 

9 a p = < 6*(d a 4>), dpcj) > 

= <d a <t> + tV da(t> N 7 d f 3<t>> 

= d a (f> ■ dfj<p + t d a n ■ dp<j) 

(1) = h aP + tb al3 . 

An important observation which will be useful later is: 

(2) VxN = V e ,xN 

for any tangent vector X. Indeed, let a = 6 o a. Then 9*X = a' and 

lim N{ ° (e)) N{ ° m 
lim N{a{e)) - N(am 

e^O e 

= VxN 

We denote the metric of M in the chart $ by g. Its components, gj ■ = di<& ■ dj$, 
are 

ga/3 = ®a a Q(j(} 
g 3 o = 
§33 = 1 

The collection of parallel surfaces St forms a foliation of M. This foliation carries 
the induced connection V given by 

V X Y = V X Y- < V X Y, N > N 
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for any two vector fields X and Y which arc tangent to the foliation. The second 
fundamental form of this foliation is the symmetric bilinear form acting on tangent 
vectors defined by 

B(X,Y) =<V X N,Y>. 

Thus 

V X Y = V X Y + B(X,Y)N. 
The components of B in the chart <£> are 

B a(i = ($>*B)(d a , dp) 

= B(d a *,df,Q) 

= <V da ^N,d/3^> 

= < Ve,(d a 4,)N,df3<f> > 

= < V da <t>N,df)<& > 

= d a n-{9ffd a 4>) 

(3) = e a a b a!i . 

The map 9 plays an important role because it allows us to extend any tensor A on 
the base surface S to a tensor 9* A on the whole bundle M. 

We conclude the description of the geometry of M with its volume element dv. 
It can be decomposed as follows 

dv = dA t dt, 

where dA t is the area element of St- Since St = 0(S), we have 

dA t = det(6»*) dA . 

Notice than in R 3 

det(0») = det(h + tb) 

= l + (trb)t+(dctb)t 2 
= l + Ht + Kt 2 , 

where H and K denote the mean curvature and the Gaussian curvature of the 
surface S respectively. 

2.2. The Infinitesimal Deformation of an Elastic Bundle. We shall assume 
that the deformation field of the base surface S is given by 

V> = ujN + W, 

where W — W <J d a <t> is an arbitrary vector field tangent to S and uu is an arbitrary 
function on S. We extend u> to the whole bundle M by the normal flow, but we'll 
continue to write uj instead of O^u — u o 9. In this section we will show that the 
corresponding deformation vector field of the elastic bundle M is given by 

(4) V = wN + 6.W-tT( u \ 
where 

T^( qi ,q 2 ,t)=T^( qi ,q 2 ) = h^d^d^. 
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From the assumption that the deformations should preserve the bundle structure 
it follows that <£> e has the following form: 



+ tn e , 



where cb e is a coordinate system for the deformed base space and n e is the normal 
to it. We have cb — cj>, n = n and 



V 



d 
de 



e=0 



(5) 

On one hand 
(6) 

while on the other hand 
d 



de 



e=0 







1 


d 


(- 


n e 1 


• n = — 




\de 


e=0 / 


2 


1c 



(n e -n e ) = 0, 



e=0 



dc 



€=0 



n t ■ d a c, 



d_ 

de 



(n e • d a cb e ) - n ■ — 



e=0 



e=0 

= — n ■ d a ip 

= —d a (ip ■ n) + ijj ■ d a n 

= -d a u + ip ■ ba^daCb 

(7) = -d a uj + b a a W a . 

Using JHJ) and Q in JSJ) we have 

V = LuN + W + ti-daUJ + b^W^h^dr^ 

= luN + {h a a + tb a a ) W a h aT d T (j) - t d a uj h aT d T cb 
= ujN + 9 a a W a h aT d T ab-td a uh aT d T cb 
= luN + 6*W -tT { ^ , 

which is what was claimed in J3J. 

2.3. The Strain Tensor. Let ~g = S\m be the metric of the bundle M. The strain 
tensor corresponding to the deformation vector field V is the symmetric tensor e 
defined by 

the Lie derivative of g in the direction of V. We now show that for any tangent 
vector field X, 

(8) e(X,N)=0. 

Let g e = $*(g) be the metric of the deformed bundle in the chart $ £ . Its components 
are: 

(g e )ij = d^e ■ dj$ e . 
The components of the strain tensor in the chart $ are 

1 d 



2 de 



e=0 
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We have 



e33 = <9 3 <I> 



de 



de 



c=0 



(9) 
and 



= 



e3o 



1 d 

2 de 



(ra e ■ d a 4> t + tn e • d a n e ) 



-=o 



(10) 



= -{-{driJ ■ n)h T ' y d a 4) ■ d a <p + n ■ d a ^p) 
= 0, 



which shows (JSJ. We can now express the strain tensor in terms of the function u> 
and the vector field W: 

1 r - 

e = 2 v g 

— 2 L -ujN+9,W-tT("> 9 

Here A stands for the symmetrization of A: 

Afj,v — —(A^u + A y ^). 
For tangential fields X, Y we have 

C aN g(X,Y)=uN < X,Y > - < C uN X,Y > - < X,C uN Y > . 

Since 

C^nX = wTVX - X(wN) 
= luC n X - (Xlu)N 

we have 

< C^nX, Y >= lo < C N X, Y > 

and 

C uN g(X,Y)=LjC N g{X,Y). 
Finally the strain tensor is given by 

(11) e = wB+v]^)-tvTH. 
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For future reference we obtain the expressions for the components of the tensors 
VT^) and V(0*W) in the chart $: 

= e a °e p T < V 9a4 ,T^\ d T <$>> 

= 9/W w -M 

(12) = 6 a a 9 fj T \7 a V T Lo, 

where V CT V T w are the components of the Hessian of to on S in the chart (j), and 

v(e*w) a p = <v a M&*w) 1 d p ^> 

= e a °e (i T < Vg^w), d T <j> > 
= e a °9 p T < v a ^(W0/^), d T <j> > 

(13) = a a 9p T T »iy a W») + 8 ct a 0p T (V a 9» T )W' 1 
Substituting lfT2^l and (fO| in (fTTJl we obtain 

(14) e aP = B af3 Lu - t9 a a 9p T W a V T ^ + 9 a °g f / (V * a W T ) + B a a 6 p T {V^) W„ . 

2.4. The Plane Stress Assumption. In linear elasticity the stress tensor a is 
related to the strain tensor e by the generalized Hooke's law: 

a ij = C ijki eH 

where 

Cijki = Cijki(qi,q2,t) 
denote the components of the elasticity 4- tensor C in the chart For a homoge- 
neous and isotropic material the elasticity tensor is of the form 

(is) c llkl = x g ij g kl + fj, g l V 7 + » s ll s Jk ■ 

In this case Hooke's law of linear elasticity takes the form 

(16) (Tij = X g mn e mn g^ + 2fi dj = X tr(e) g tj + 2fi e l3 . 
The strain energy corresponding to the strain e is defined by 

£ = \ f C(e,e)dv=l f C^ kl e l3 e kl 

It is now possible to derive a shell theory based solely on the Kirchhoff-Love hy- 
pothesis using the expression for the strain tensor l|14|) . It turns out however that 
a more realistic shell model is obtained using the assumption of plane stress: 

(17) <7 33 = 0. 

In linear elasticity this assumption contradicts the Kirchhoff-Love hypothesis. We 
follow a standard approach in modifying the Kirchhoff-Love hypothesis (see jjj): 
we no longer assume that 

e33 = 0, 

but we continue to assume that 

e 3a = 0. 
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Let e denote the tangential part of the strain e. From (|16f) and l|17f) it follows that 
= 0-33 = Ag m "e m „ + 2/xe 33 = Ag"%„ + (A + 2/x)e 3 3 

and therefore 

e33 = -AT2^ g " e - = -AT^^- 
We use the last equation to simplify the expression for the strain energy: 

C(e, e) = A (g mn e mn ) 2 + 2fig mn g kj e mj e kn = A(tre) 2 + 2 A/ tr(e 2 ) . 

Since 

tee = tre + e 33 = (1 - j^) tre = _ tre 



and 



tr(e 2 ) = tr(e 2 ) + e 2 3 = tr(e 2 ) + - (tre) 2 



we obtain 

C7(e,e) = 4A I ^(tr,) 2 + 2, Ix ^ F (tr.) 2 + 2 M tr(. 2 ) 
" A/i -(tre) 2 + 2//tr(e 2 ) . 



A + 2/x 

2.5. Variation of the Strain Energy. In order to obtain the equilibrium equa- 
tions we proceed to calculate the variation of the strain energy. Let uj be an arbitrary 
function and W an arbitrary vector field on S. The corresponding variation of 
the energy is 

.. £(w + eu>, W + eW) - £(w, W) 
ob = hm . 

Let e = e(V) be the strain corresponding to the deformation V. Since 

e(V + V) = e(V) + e(V) 
the variation of the strain energy is 

X + 2/j, J M J M 



h ° ' A^ 



mJ-h \ A + 2 M 

= f [ ° A a ^ s e af3 e jS dtdA , 
Jn J -h 

where we define 

A a ^ s ( qu q 2 ,t)= L^—g^ + ^g^g^ det(^). 

For the calculations that follow it is useful to note the following symmetries of A: 
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We use ()14JI and integrate by parts omitting the boundary terms to rewrite the 
energy variation as follows 



JQ J -ho V ' 

= f f ° {K^ s e a0 B l5 -tV T V a {k al3 '' s e a ^ a e 5 T )) CodtdA , 
Jn J-h 

+ 11° (A af * e e af} 6S0 s T V tr 6 T ' t - V CT (A^^A's/)) W„ dtdA Q , 
Jn J -h 



Collecting the terms we obtain the equations for the normal and the tangential 
components of the force 



/ 3 = [ ° {A a ^ 5 e a pB l5 - t V r V CT (A a ^ 5 e a ^6 s T )) dt 

J -ho 

f = / ° {A aPlS e a0 6^6 & T V a e^ - V ff (A^e^/g/)) dt 

J -h a 



Using 114(1 again we can finally bring these equations into the following form 



(18) / 3 = Auj + y a y T {A M V m V„cj) - V a V T {A uj)-A V a V T Lu 

+ ® v w v + $ flv v„w v - V ff V r (^ CTT ^) - v M v„(* ffT/w v« 7 w r ) 

(19) /" = W W v + ?T™ V a W T - Wa^f^Wr) - V^Q^VaWr) 
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where the coefficients are defined as follows: 

(20) A = J A a ^ s B a( 3B lS dt 

(21) T T ^ = J k al3 '< 5 e a Bf } T e^e & l 't 2 dt 

(22) = J A a ^ 5 B a ^e s u tdt 

(23) ^ = J K^B^e^OsPiVre^dt 

= A V T 6/ 

(24) ^ = J k a0 ~> s B aP e^e s T 9 T v dt 

(25) = J h al3 '< s e a ' J e p T {v ^^e^Os" t 



dt 



A^ T ^V a b T » 



(26) tf"""" = j A^/g/fl/Virft 

(27) = J A a i 3 "' s e ol a 6p T {y (7 e T ' J -)e 1 x 6 s p {Vx0 P v )dt 

= A° TXp V a b T »V x b p » 

(28) VT" = J A a ^ s 6 a Pg p v 6 1 T 6 s x (V T 6 x P)dt 

= * V„6/ 

(29) If™" = J K a ^ s e^ g0 T e^g s v dt 

2.6. The Plate Equation. For a flat shell surface S the only non-zero coefficients 
are A and £1 is zero because it is an integral of an odd function on a symmetric 
interval), and 

A ot/» = A/i h^h^ + uh' Tfi h TV 
3 

n - 2/to A 0-1 ""*' 
The equations {TBI - CHJ reduce to 

/norma/ 

f tangential = Ih^D \7 ~ W 

where 

D = 2/j, 



2 2 

/normal = 7-7^0^^ ^ 



A + 2/Z 

The first equation is the classical equation of an elastic plate (see The second 
equation does not appear in classical plate analysis where it is assumed that the 
plate undergoes only vertical motion and the shearing forces are negligible. 
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3. The Immersed Boundary Method 

The immersed boundary method of Peskin and McQueen is a general numerical 
method for modeling an elastic material immersed in a viscous incompressible fluid. 
Typically the immersed material has been modeled as a collection of fibers. For 
details and references to many applications see ^3] and This section outlines 
the immersed boundary method with the elastic material being modeled as a shell. 

3.1. The Equations of the Model. For the description of the fluid we adopt a 
standard cartesian coordinate system on R 3 . The immersed material is described 
in a different curvilinear coordinate system. The equations can be partitioned into 
three groups: the Navier Stokes equations of a viscous incompressible fluid, the 
equations describing the elastic material and the interaction equations. Accordingly 
there are two distinct computational grids: one for the fluid and the other for the 
immersed material. The purpose of the interaction equations is to communicate 
between these two grids. 

We first turn to the Navier-Stokes equations. Let p and fj, denote the density 
and the viscosity of the fluid, and let u(x, t) and p(x, t) denote its velocity and 
pressure, respectively. Then 

(30) p^+u-Vu^j = -Vp + fiV 2 u + F 

(31) Vu = 

where F denotes the density of the body force acting on the fluid. For example, if 
the immersed material is modeled as a thin shell, then F is a singular vector field, 
which is zero everywhere, except possibly on the surface representing the shell. 

Let X(q, t) denote the position of the immersed material in R 3 . For a shell, q 
takes values in a domain C R 2 , and X(q, t) is a 1-parameter family of surfaces 
indexed by t, i.e., X(q, t) is the middle surface of the shell at time t. Let f(q, t) 
denote the force density that the immersed material applies on the fluid. Then 

(32) F(x, t) = j f (q, i)<5(x - X(q, <)) dq, 

where <5 is the Dirac delta function on R 3 . This equation merely says that the fluid 
feels the force that the material exerts on it, but it is important in the numerical 
method where it is one of the interaction equations mentioned above. The other 
interaction equation is the no slip condition for a viscous fluid: 

— = u(X(q,*),i) 

(33) = J u(x,i)<5(x-X(q,<))dx 

We will now complete the description of the system by writing down the third 
group of equations, the equations describing the shell. Let Xo(q) denote the initial 
position, which we take as our equilibrium reference configuration, of the middle 
surface of the shell, and let T M = T M (q) = g|- (/i = 1, 2) be its tangent coordinate 
vector fields. The displacement from the stationary configuration can be described 
in terms of the function w(q, t) and the vector field W(q, t) defined by the following 
equation: 

X(q, t) - X (q) = w(q, t)N(q) + W(q, t), 
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where N = N(q) denotes the normal to the surface Xo(q) and W = W^T^ is 
tangent to Xo. We decompose the force f into its normal and tangential components 
as well: 

f = / 3 N + /"T„ 

The components /* are related to u> and by the shell equations (|18|l (|19ll . 
This completes the description of the fluid-shell system. 

3.2. The Numerical Method. Let At denote the duration of a time step. It will 
be convenient to denote the time step by the superscript. For example u n (x) = 
u(x, nAt). At the beginning of the n-th time step X™ and u™ are known. Each 
time step proceeds as follows. 

(1) Compute the force f" that the shell applies to the fluid. 

(2) Use ().'>2t to compute the external force on the fluid F". 

(3) Compute the new fluid velocity u ra+1 from the Navier Stokes equations. 

(4) Use (|33f) to compute the new position X ra+1 of the immersed material. 
The computation of the force in step 1 is explained in detail in the next section. 
We shall now describe in detail the computations in steps 2 — 4, beginning with 
the Navier-Stokes equations. 

The fluid equations are discretized on a rectangular lattice of mesh width h. We 
will make use of the following difference operators which act on functions defined 
on this lattice: 

, </f>(x + hei) - 0(x) 



h 

(x) - 0(x - he { ) 
h 

(x + hei) ~ 0( x — hei 



(34) D+0(x) = 

(35) D7>(x) = 

(36) £>?0(x) = , )/; 

(37) D° = (DlDlDl) 

where i = 1, 2, 3, and ei, e2, e3 form an orthonormal basis of R 3 . 

In step 3 we use the already known u™ and F" to compute u" +1 and p n+1 by 
solving the following linear system of equations: 

(n+l n ^ \ 3 

fc=i / fc=i 

(39) D°-u" +1 = 

Here u^D^ stands for upwind differencing: 

± _ f ulD- k ut > 

U k^k ~ I u n D + u n < 

Equations (|38|l - (|39|) are linear constant coefficient difference equations and, there- 
fore, can be solved efficiently with the use of the Fast Fourier Transform algorithm. 

We now turn to the discretization of equations (|32|l . (|33|l . Let us assume, for 
simplicity, that VI C R 2 is a rectangular domain over which all of the quantities 
related to the shell are defined. We will assume that this domain is discretized with 
mesh widths Aqi , Aq2 and the computational lattice for £1 is the set 

Q = {(iiAqi,i 2 Aq 2 ) | h = 1 . . .n\, h = 1 . . . n 2 } . 



MODELING ELASTIC SHELLS IMMERSED IN FLUID 



15 



In step 2 the force F" is computed using the following equation. 

(40) F"(x) = f"(q)^(x - X"(q))Aq 

qeQ 

where Aq = Aqi Aq2 and Sh is a smoothed approximation to the Dirac delta func- 
tion on R 3 described below. 

Similarly, in step 4 updating the position of the immersed material X n+1 is done 
using the equation 

(41) X" +1 (q) =X"(q) + Ai^u"+ 1 (x)^(x-X"(q))/i 3 , 

X 

where the summation is over the lattice x = (hi, hj, hk), where i, j and k are 
integers. 

The function Sh which is used in l(4()jl and 141(1 . is defined as follows: 
where 

( |(3-2|r| + Vl + 4|r| -4r 2 ) M < 1 
4>(r) = l i-0(2-|r|) l<|r|<2 

[ 2<|r| 

For an explanation of the construction of Sh see |13) . 

3.3. Computation of the Elastic Force. The force f™, that the shell applies to 
the fluid is computed by discretizing equations (|18fl - (|19|l . We shall now describe 
how to discretize these equations and how to compute various geometric quantities, 
such as the Christoffel symbols, the second fundamental form, etc. These basic 
geometric quantities, as well as the coefficients (|20|l - (|29f) that depend on them, 
can be computed once in the initialization step of the algorithm and stored for 
subsequent use. 

A covariant derivative of a (p, g)-tensor A is defined by 

q p 

y<±6) V Q ^ Ml ... Mp — u aS± p , 1 ...p lp T 1 acr^Vi...^, / , L an„, Hi ...<J...U.„ ' 

fc=l m=l 

where 

^ \iv = 2 9 (9^<y,v "I" 9ov-V- ~ 9fJ.v-Xr) 

are the Christoffel symbols, t/ M „ = T„ • T„ is the metric of the surface Xo, is 
the inverse of <7 M „, and comma denotes partial differentiation (i.e., jQ = d a <f). To 
simplify the discretization of the force equations we introduce a single difference 
operator D a which uses center differencing in the interior of the domain and either 
forward or backward differencing on the boundary: 

( D+(j)(q) q a = Aq a 
D a <t){q) = 1 £>°0(q) 2Aq a < q a < (n a - l)Aq a 
{ D~4>(c\) q a = n a Aq a 

Here the operators D+ , Z}° and D~ are defined on the lattice Q analogously to the 
definitions IpTH) - (|3TJl . We define the discrete covariant derivative of a (p 7 g)-tensor 
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A by 

q p 
k—1 m— 1 

The following computations are carried out in the initialization stage of the 
algorithm: 

(1) Compute the tangent vector fields T a = D a Xo. 

(2) Compute the unit normal vector field N = ^ XT2I ' 

(3) Compute the metric g^ v — T M • T„ and its inverse g^" . 

(4) Compute the second fundamental form: b^ v = Z? P N ■ T„ . 

(5) Compute the Christoffel symbols: 

r£„ = \g aX {D v g^ + D^g av - D a9fw ) . 

(6) Compute the derivative of the second fundamental form: 

B a bff = D a b^ + Tl a bf ~ T^b^ . 

(7) Compute the coefficients (1201) - (J29JI . I n practice instead of evaluating the 
integrals in these expressions it is easier to assume that the shell is suffi- 
ciently thin and to drop the terms involving high powers of the thickness. 

The computation of the force during each time step proceeds as follows: 

(1) Decompose the displacement into its tangential and normal components: 

lo = (X" - X ) • N 

W a = (X" - Xo) • T Q 

(2) Compute D^lo, the components of du), the 1-form dual to the gradient of 

LO. 

(3) Compute D^D^lo, the components of the hessian of to. 

(4) Compute the components D^W V . 

(5) Compute the normal and the tangential components of the force using the 
following discretization of the equations Ijl8|l - l|19|) : 

f = Alo + D a D T {A D^DyLo) — D a D T (A lo) - A D a D T LO 

(44) + <$> V W V +^ U D^W, - D a D T (^ aT W^) - DM^ aTtlu D a W T ) 

= W v W v + ff ™ D a W T - D„(fF" T W T ) - Dviff^DtrWr) 

(45) + $^ lo - D U (^ V " lo) - y^ T D a D T uj + b v {^ ViiaT D a D T uo) 

To complete the computation, we express the force in cartesian coordinates using 

f=/ 3 N + /%. 

4. An Application: Modeling a Shell Immersed in Fluid 

In this chapter we present a model shell which is a prototype of a piece of the 
basilar membrane of the cochlea. The immersed boundary method for shells was 
implemented using the C programming language and numerical experiments were 
carried out with this shell. We describe the model that was used, examine the 
convergence of the algorithm in this case and discuss the results. 
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n — 1 cm 


IfiT-ip-fV) of trio sirlo of flip flnirl rnV)p 


AT OO CiA 1 9Q 

iv — oz, 04, izo 


size of the fluid grid is 


h — a 


fluid mesh width 


p = 1.034 cm g- 3 


11U1Q (.ICIlolLy 


v = 0.0197 g cm^ 1 sec -1 


fluid viscosity 


Lbm = 3.5 cm 


length of the basilar membrane 


L = 0.5 cm 


length of the model shell 


wq — 0.015 cm 


width of the basilar membrane at the base 


w\ — 0.056 cm 


width of the basilar membrane at the apex 


w{qx) =w + L ^ M (wi w ) 


width of the model shell 


%i) w 0.001(1.0+5.0 gi) cm 


thickness of the model shell 


a = 1.8-ir/L 
R = cm 






i? = 0.01 cm 




ni = 1280^ 


first dimension of the surface grid 


n 2 = 48^ 
A ?i - ni -i 


second dimension of the surface grid 


surface mesh width 




surface mesh width 


A = 26197503.0 g cm" 1 sec~ 2 


first Lame coefficient 


fi = 523950.0 g cm- 1 sec" 2 


second Lame coefficient 


At = 0.5, 1.0, 2.0, 4.0 x 10" 8 sec 


time step 


T = 2.0 x 10 -6 sec 


total simulation time 



4.1. The Model Shell. We construct a shell similar to a piece of the basilar 
membrane. Our shell will be approximately one seventh of the length of the basilar 
membrane in the human cochlea and it will have more curvature. Introducing 
more curvature allows us to pack a longer shell into a cube of fluid of a given size. 
Choosing a small fluid cube enables us to achieve a better numerical resolution of 
the fluid. 

The surface is a narrow helicoidal strip defined as follows. Consider the curve 

7(i) = (i?cos(ai), R sin(crf), Hat) 

where R, a and H are constants. They are specified, along with the other param- 
eters that we use below, in Tabled The vector field 

N(t) = (- cos(ai), - sin(crf), 0) 

is a unit normal field along the curve 7. We define the surface of the shell by the 
following equation: 

X(gi,q 2 ) - 7(31) + (32 - ^f^j N( Q1 ), 

< qi < L 7 < q 2 < — - — 
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0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



Figure 1 . The thickness of the model shell. 

The width w is chosen to grow linearly, similar to the width of the basilar membrane 
(see Table 0). This surface is discretized on a grid of size ri\ x n 2 as follows: 

X fel , fc2 = 7 (fc 1 A gl ) + (fe 2 Ag 2 (fc!A gl ) - w ^ q ^ AT( fcl A gl ) 

fci = 1, ...,m, k 2 = 1, n 2 . 

The compliance of an elastic shell is defined as the amount of volume displaced 
per unit pressure difference across the shell. Von Bekesy found that the compli- 
ance per unit length of the basilar membrane varies exponentially along the basilar 
membrane as e cqi , where c -1 = 0.7 cm (see To achieve this compliance in the 

model shell we use the shell equations to estimate the required thickness. (Notice 
that nothing changes in the derivation of the shell equations when the thickness ho 
is assumed to be a function on the middle surface, rather than being a constant. 
The thickness enters the equations only in the definition of the coefficients 1201) - 
(|29Jl ). We obtain an estimated thickness 

h(qx) = 0.001 ^1 + ~qi\ " 1(H ?1 cm, 

which on the interval < qi < 0.5 is very close to a straight line (Figure QJ. The 
above choices of the width, the thickness and the Lame coefficients (see Table 
yield an estimated compliance similar to the one measured by von Bekesy for the 
basilar membrane. 

4.2. The Numerical Experiments. The immersed boundary method for shells 
was implemented using the C programming language. The program ran on the 
Cray C-90 at the Pittsburgh Supercomputing Center as well as on Silicon Graphics 
workstations. 

The numerical experiments are set up as follows. The shell X is placed in a 
periodic cube of fluid, i.e., a three-dimensional torus. The length of this cube's side 
is a = 0.1 cm, so it is much smaller than the cochlea, whose volume is approximately 
1 cm 3 . The fluid density p, and the fluid viscosity v, are chosen equal to the 
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corresponding parameters of the cochlear fluid. The fluid equations are discretized 
on a cubic grid of N 3 points and the grid size for the shell is chosen correspondingly 
such that its mesh widths, Aqi and Aq2, are approximately equal to half of the 
fluid mesh width h (Table 0) . This choice is necessary to prevent the fluid from 
leaking through the shell. The edges of the shell's boundary are clamped to a fixed 
location in space with two rows of springs along each edge. 

At time the system is perturbed with a force impulse in the fluid and the simu- 
lation is run for a fixed period of time Tq = 2.0 x 1CP 6 sec. The force impulse acting 
at time is represented by a constant singular vector field in the vertical direction 
defined on the horizontal plane z = 0.1 cm (where z is the vertical coordinate) 
with the force density of 4 x 10~ 7 g cm _1 scc -2 . Since the force source is broadly 
distributed through a horizontal plane within the fluid, wave propagation within 
the basilar membrane can only arise as a consequence of its material properties and 
cannot be the result of the location of the force source. 

Immersed boundary computations typically require large scale computing re- 
sources. Because of time and storage limitations, it is not possible to conduct 
an extensive empirical study of the algorithm's convergence. At present it is not 
practical to implement the method with a fluid grid of more than 128 3 points. 
The experiment has been repeated with N = 32, 64 and 128 and time steps At = 
0.5, 1.0, 2.0 and 4.0 x 10 -8 sec. Numerical stability condition forces a choice of 
such small time steps. On the other hand, reducing the time step further may result 
in a significant machine precision error. 



4.3. Convergence of the Algorithm. Let Xi(t) and X 2 (t) denote the position 
of the shell at time t obtained in two different numerical experiments. To measure 
the relative difference at time t we calculate 



(46) 



E{t) 



|X 1 (t)-X 2 (t)|, 



IXiW-X^lp 
where the norms are L p -norms with p = 1, 2 and oo. 



Relative error as a function of time step 



Top graph: relative max-norm error 
Middle graph: relative L2 norm error 
Bottom graph: relative L1-norm error 




Figure 2. The relative difference between Xx,i28 an d X 2j i28- 
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For any two experiments the function E has been found to decrease initially very 
sharply, stabilizing in the second half of the time interval. Two typical graphs are 
shown in Figure El and Figure |3 Therefore, we will measure the difference between 
two computed solutions on the time interval [-y-, To] using the following space-time 
norms: 

||X 1 -X 2 || P = Yl |Xi(t)-X 2 (<)| p 

te[0.5 T , To] 

where the summation is taken over the set of 25 time values common to all of the 
tests performed and the L p -norms are computed on the "common" grid of 257 x 12 
points. 



Relative error as a function of time step 




50 100 150 200 250 300 350 400 

time step 



Figure 3. The relative difference between Xo. 5,12s and Xo.5,64- 

The computations were performed with double precision floating point on a 
Silicon Graphics computer. The results of the measurements are summarized in 
Table El and Table 01 These results indicate that the algorithm converges when the 
time step At is linearly proportional to the mesh width h and both tend to zero. 
Let X 5 denote the position of the shell computed with At = S x 10 -8 and N = N. 

Table 2. Norm comparison between different numerical tests. 



Tests compared 


L 1 -norm 


L 2 -norm 


L 00 -norm 


1.0/128 2.0/64 
2.0/64 4.0/32 

0.5/128 1.0/64 
1.0/64 2.0/32 


9.6290 x 10~ U8 
2.3984 x 10~ 07 
9.2989 x 10~ 08 
2.4082 x 10~ 07 


1.8547 x 10~ U ' J 
4.6353 x 10~ 09 
1.8086 x 10~ 09 
4.7011 x 10-° 9 


9.7749 x 10" 11 
1.8326 x 10- 10 
9.5158 x lO" 11 
1.9287 x 10~ 10 



Using the data in Table El we have two estimates of the order of convergence of the 
algorithm: 

r _ 1 Q ( l|Xl,128 — X 2! 64||p \ 
V ||X2,64 — X4 i 3 2 || p / 
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Table 3. Convergence rate estimates. 





i 1 -norm 


L 2 -norm 


I/°°-norm 


T2 


1.3729 
1.3166 


1.3781 
1.3215 


1.0192 
0.9068 



and 

r _ j Q ( I|Xq.5,128 ~ Xi,64||p | 
2 \ ||Xl,64 — X 2 ,32||p / 

The values of n and r 2 are given in Table 01 As the time step becomes smaller, 
the machine precision error becomes more significant leading to a slower rate of 
convergence. This is apparently the reason that we have r 2 < r% . 

4.4. The Traveling Wave. The test model described above was run with N — 
128 and time step At — 1.0 x 10 -8 seconds for a total of more than 2400 time 
steps. Although this model is still far from a realistic model of the cochlea, several 
qualitative features, characteristic of the cochlear wave mechanics, were already 
observed in this experiment. 

In the experiment a traveling wave is produced in response to the force impulse 
in the fluid. The wave propagates in the direction of increasing compliance within 
the elastic shell. This agrees with von Bekesy's observation that the traveling wave 
in the cochlea always propagates from the base to the apex (see JHl)- Snapshots 
of this wave are shown in Figure The ten snapshots were taken every 200 time 
steps beginning with time step 600. The first five appear in the first column, the 
rest in the second. Since the displacements of the shell are too small to be seen 
with the naked eye, they are represented in the pictures on a gray scale. Black 
color indicates the maximal possible displacement down, and white, the maximal 
displacement up (the scale is symmetric with respect to zero). The initial force 
impulse has pushed the fluid down and, as can be seen from the first snapshot, 
after 600 time steps the shell is displaced downwards. In the following snapshots 
we can see that the stiffer part of the shell near the base is bouncing back and a 
wave starts propagating towards the apex. 

Conclusion and Further Research 

This paper describes an extension of the immersed boundary method to elastic 
shells in a viscous incompressible fluid. The method was developed as a part of a 
project to construct a three-dimensional computational model of the cochlea. It 
is based on shell equations derived using techniques of differential geometry. The 
resulting method is a practical method which has been implemented and tested on 
a prototype of a piece of the basilar membrane. We have examined the convergence 
of the algorithm in this case. Numerical experiments indicate that the algorithm 
has the first order convergence rate when the time step is chosen to be linearly 
proportional to the fluid mesh width. 

The numerical experiments have shown a traveling wave propagating in the test 
model shell in the direction of increasing material compliance in response to external 
impulsive excitation of the fluid. This reproduces the so-called "travelling wave 
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Figure 4. The traveling wave in the test shell. 

paradox" observed in the cochlea: the traveling wave always travels from the base 
to the apex of the cochlea independently of the location of the impulse source. 

The numerical method described here was subsequently used in a construction 
of a complete three-dimensional computational model of the macro-mechanics of 
the cochlea. Additional difficulties in the construction of this model are caused 
by the large scale of the immersed boundary computations required, and by the 
presence in the fluid of several different elastic materials in addition to the basilar 
membrane. The size of a human cochlea is on the order of 1 cm x 1 cm x 1 cm, 
while the basilar membrane is about 3.6 cm long and is very narrow (150 — 560 
/an). Therefore, larger fluid grid is necessary in the full cochlea simulation in order 
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to resolve the 100 /j,m scale corresponding to the width of the basilar membrane. 
Immersed boundary computations require large scale computing resources: a fluid 
grid of 256 x 256 x 128 points needed for the full cochlea model requires a significant 
amount of computer memory. The CFL condition and the stability conditions 
imposed by the stiffness forces of the immersed boundaries force a choice of a very 
small time step (approximately 50 nano-seconds) when the fluid mesh width is 
small. The convergence study carried out in this paper indicates that decreasing 
the fluid mesh width by a factor of two necessitates a corresponding decrease in the 
time step by approximately a factor of two. This has indeed been further verified 
in the construction of the complete cochlea model. 

The construction of the cochlea model has been achieved in collaboration with 
Julian Bunn on Hewlett-Packard computers at the Center for Advanced Compu- 
tational Research at Caltech. The results of this work will be described in future 
publications. 
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